function VC_IV()
% Calculate IVs based on data passed through DFILE, CONTROL(sf), current and voltage
% plot_flag : 0 is to display only,1 is to print, and 2 is to print and then close
% Data is assumed to have been read and corrected when it is passed to this routine
% Step times for for the data are derived from contents of the CONTROL matrix
% rather than the DFILE values, so we can correct them if necessary.
testflag = 1;

QueMessage('VC analysis', 1); % clear the que
htag = findobj('Tag', 'CTL');
sf = get(htag, 'Value');
if(sf <= 0) return;
end
pflag = getplotflag;

if(testflag)
    for i = 1:length(sf)
        vc_iv2(sf(i), pflag); % handles iv, iva, ivp and the like.
    end;
    return;
else

    try
        for i = 1:length(sf)
            vc_iv2(sf(i), pflag); % handles iv, iva, ivp and the like.
        end;
    catch
        watchoff;
        QueMessage('Error in VC Analysis routine', 1);
    end;
    return;
end;
   

function vc_iv2(sf, plot_flag)

global VOLTAGE CURRENT DFILE
global CONTROL

ctl = 0;

[DFILE] = synch_file(DFILE, sf);
if(isempty(DFILE) | isempty(VOLTAGE))
   fprintf('? could not get dfile');
   return;
end;

protocol=lower(CONTROL(sf).protocol);
QueMessage(sprintf('Processing %s, %s for IV', DFILE.filename, protocol));
RATES = DFILE.rate.*DFILE.nr_channel/1000;

[records,pts]=size(CURRENT);
time=make_time(DFILE);
tmax=max(max(time'));

% mean times for standard IV are computed as percentages from the #2 window
ts1=number_arg(CONTROL(sf).durho);
ts2=number_arg(CONTROL(sf).durs1);
ts3=number_arg(CONTROL(sf).durs2);
ts4=number_arg(CONTROL(sf).durs3);
tdeadwin = number_arg(CONTROL(sf).deadwin);
ACT=1;
INACT=2;
BOTH=3;
analysis=0; % none set.
tailanalysis = 0;
protocol=lower(deblank(strjust(protocol, 'left'))); % remove leading and trailing spaces from the string.

% set up parameters for analysis
switch(protocol)
case {'iv', 'vciv', 'pativ', 'iv2', 'vciv2', 'vc_iv2'} % standard IV from -60 mV
   analysis=ACT;
    tailanalysis = 1;
    tdelA=ts1; tstepdurA=ts2;
   
case 'iva' % IV obtained after a prepulse step
   analysis=ACT;
   tailanalysis = 1;
   tdelA=ts1+ts2; tstepdurA=ts3;
case {'inact', 'ivp', 'mc-inact'} % inactivation of I after varying V during prepulse step
   analysis=INACT;
   tdelI=ts1+ts2; tstepdurI=ts3;
   tv0=ts1+ts2*0.5;
   tv1=ts1+ts2*0.9; % voltage measured at a different time!
case {'qinactp2','qinactp1','mc-qinactp2','mc-qinactp3'} % inactivation of I on step 3 after varying V during step 2
   analysis=BOTH;
   tdelA=ts1+ts2+tdeadwin; tstepdurA=ts3; % set activation window
   tdelI=ts1+ts2+ts3+tdeadwin; tstepdurI=ts4/2; % set inactivation window
   tv0=ts1+ts2+ts3*0.5;
   tv1=ts1+ts2+ts3*0.9; % voltage measured at a different time!
case 'qinactp3' % inactivation of I on step 3 after varying V during step 2
   analysis=BOTH;
   tdelA=ts1+ts2; tstepdurA=ts3; % set activation window
   tdelI=ts1+ts2+ts3; tstepdurI=ts4/2; % set inactivation window
   tv0=ts1+ts2+ts3*0.5;
   tv1=ts1+ts2+ts3*0.9; % voltage measured at a different time!
otherwise
   QueMessage(sprintf('VC_IV: processing for protocol %s not implemented', upper(protocol)));
   return;
end

% compute windows for activation analysis
if(bitget(analysis,1)) % any case where activation is done...
   m1A=tdelA+tstepdurA*0.95;
   m2A=tdelA+tstepdurA*0.99;
   %peak times:
   p1A=tdelA+tstepdurA*0.005;
   p2A=tdelA+tstepdurA*0.5;
   
   mean_t1A=floor(m1A./RATES);
   mean_t2A=floor(m2A./RATES);
   peak_t1A=floor(p1A./RATES);
   peak_t2A=floor(p2A./RATES);
end

if(bitget(analysis,2)) % prepulse voltage measure
   mean_v0=floor(tv0./RATES);
   mean_v1=floor(tv1./RATES);
   m1I=tdelI+tstepdurI*0.95;
   m2I=tdelI+tstepdurI*0.99;
   %peak times:
   p1I=tdelI+tstepdurI*0.005;
   p2I=tdelI+tstepdurI*0.5;
   mean_t1I=floor(m1I./RATES);
   mean_t2I=floor(m2I./RATES);
   peak_t1I=floor(p1I./RATES);
   peak_t2I=floor(p2I./RATES);
end

% Smooth the current out a bit
for i=1:records
   fsamp = 1000/RATES(i); % get sampling frequency
   fco = 2000;		% cutoff frequency in Hz
   wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
   if(wco < 1) % if wco is > 1 then this is not a filter!
      [b, a] = fir1(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
      current_s(i,:) = filtfilt(b, a, CURRENT(i,:)); % filter all the traces...
   end
end

% now compute the relevant IVs.
% there are up to 4 possible IVs:
% volt_ssA is the activation voltage; volt_ssI is the voltage for inactivation protocols.
% curr_peak_A is peak current for activation protocols; curr_ss_A is the peak current
% curr_peak_I is peak current fir inactivation protocols; curr_ss_I is the steady state current

   T=32.0;	% temp, deg C
   VR=-81.5;

if(bitget(analysis, 1)) % compute activation IVs, peak and ss, with current and voltage measured at the same time
   QueMessage('Computing activation');
   tic;
   for i=1:records
      volt_ss_A(i,:)=mean(VOLTAGE(i,mean_t1A(i):mean_t2A(i)));
      curr_ss_A(i,:)=mean(current_s(i,mean_t1A(i):mean_t2A(i))); % current and voltage are measured at the same time.
      [curr_peak_A(i), peak_i(i)]=max(current_s(i,peak_t1A(i):peak_t2A(i)));
      % recompute the peak current getting a time window of 1 msec around the detected peak
      pk0=floor((peak_i(i)+peak_t1A(i)-1)-0.5/RATES(i));
      pk1=floor((peak_i(i)+peak_t1A(i)-1)+0.5/RATES(i));
      if(pk0 <= 1) pk0 = 1; end; if(pk0 > size(current_s, 2)) pk0 = size(current_s, 2) -1; end;
      if(pk1 <= 1) pk1 = 1; end; if(pk1 > size(current_s, 2)) pk1 = size(current_s, 2) -1; end;
      curr_peak_A(i)=mean(current_s(i,pk0:pk1));
      volt_trans_A(i)=VOLTAGE(i,peak_t1A(i)+peak_i(i)-1);
   end;
   
   % sort the data so we can interpolate correctly
   
   [v_ssa, ivssa] = sort(volt_ss_A);
   c_ssa = curr_ss_A(ivssa);
   v_tra = volt_trans_A(ivssa);
   c_pka = curr_peak_A(ivssa);
   
   % interpolate to get correct values due to v sag...
   
%   c_ssa_i = interp1(v_ssa, c_ssa, v_tra, 'cubic'); % use polynomial fit for better working
   pfit = polyfit(v_ssa, c_ssa, 5);
   c_ssa_i = polyval(pfit, v_tra); % recompute it at the points corresponding to the transient
   c_tra = c_pka - c_ssa_i; % subtract the interpolated versions...
   
   maxiter = 100;
   m=max(gradient(c_tra./v_tra)); % difference current
   [fpk_a, chisqpk_a, niterpk_a, fit_pka, errpk_a, deppk_a] = mrqfit('boltzmann', ...
      [0 m 1 -20], v_tra, c_tra, [], [0 1 1 1], [], [], maxiter);
   
   act.timestamp = datestr(now);
   act.edit = 0; % not editable - analysis results...
   act.pk_imax=fpk_a(2); act.pk_k=fpk_a(3); act.pk_v05=fpk_a(4);
   act.pk_err=sqrt(sum((fit_pka-c_tra).^2))/length(fit_pka); act.pk_Er=VR; act.pk_T=T;
   act.pk_V=v_tra; act.pk_I=c_tra;
   
   m=max(gradient(curr_ss_A./volt_ss_A)); % difference current
   [fss_a, chisqss_a, niterss_a, fit_ssa, errss_a, depss_a] = mrqfit('boltzmann', ...
      [0 m 1 -20], v_ssa, c_ssa, [], [0 1 1 1], [], [], maxiter);
   act.ss_imax=fss_a(2); act.ss_k=fss_a(3); act.ss_v05=fss_a(4);
   act.ss_err=sqrt(sum((fit_ssa-c_ssa').^2))/length(fit_ssa); act.ss_Er=VR; act.ss_T=T;
   act.ss_V=v_ssa; act.ss_I=c_ssa;
   watchoff;
   fit_t=toc;
   QueMessage(sprintf('Activation IV calculations took %8.3f sec', fit_t));
   CONTROL(sf).act = act; % store activation results in its own structure.
   
end
% compute inactivation IVs with voltage measured prior to current...
% added double boltzman fit on peak current using Scott's dll version...
if(bitget(analysis, 2))
   max_negV = -140;
   max_posV = -60;
   QueMessage('Computing Inactivation');
   tic;
   for i = 1:records
      volt_ss_I(i,:)=mean(VOLTAGE(i,mean_v0(i):mean_v1(i))');
      curr_peak_I(i,:)=max(current_s(i,peak_t1I(i):peak_t2I(i))');
      curr_ss_I(i,:)=mean(current_s(i,mean_t1I(i):mean_t2I(i))');
   end;
   [v_ssi, ivssi]=sort(volt_ss_I); % sorted voltage... according to the activation list.
   c_pki=curr_peak_I(ivssi);
   c_ssi=curr_ss_I(ivssi);
   vok=find(v_ssi < -10 & v_ssi > max_negV); % only include data for steps negative to actual pulse
   alpha=(10^-3)*9.648*10^4/(8.3143*(273.16+T));	%1/25mV  =F/RT 
   m = max(curr_peak_I);
   inact.timestamp = datestr(now);
   inact.edit = 0; % not editable
   [fp2, chisq2, niter2, fit_pki2, err2, dep2] = mrqfit('boltzmann', [0 m/2 -0.1 max_negV m/2 -0.1 0], v_ssi(vok), c_pki(vok), [], [1 1 1 1 1 1 1], [], [], 250);
   [fp, chisq, niter, fit_pki, err, dep] = mrqfit('boltzmann', [0 m/2 -0.1 max_posV], ...
      v_ssi(vok), c_pki(vok)); % fit sorted voltage
   % double boltzman for peak current
   inact.pk2_imax=fp2(2); inact.pk2_k=fp2(3); inact.pk2_v05=fp2(4);
   inact.pk1_imax=fp2(5); inact.pk1_k=fp2(6); inact.pk1_v05=fp2(7);
   inact.pk2_err=sqrt(sum((fit_pki2-c_pki(vok)').^2))/length(fit_pki2); inact.pk2_Er=-81.5; inact.pk2_t=T;
   % single boltzman for peak current
   inact.pk0_imax=fp(2); inact.pk0_k=fp(3); inact.pk0_v05=fp(4);
   inact.pk0_err=sqrt(sum((fit_pki-c_pki(vok)').^2))/length(fit_pki); inact.pk0_Er=-81.5; inact.pk0_t=T;
   inact.pk_V=v_ssi; inact.pk_I=c_pki;
   
   m = max(c_ssi);
   
   % Fit the "steady-state" curents 
   [fs2, chisq2, niter2, fit_ssi2, err2, dep2] = mrqfit('boltzmann', [0 m/2 -0.1 max_negV m/2 -0.1 0], v_ssi(vok), c_ssi(vok), [], [1 1 1 1 1 1 1], [], [], 250);
   [fs, chisq, niter, fit_ssi, err, dep] = mrqfit('boltzmann', [0 m -0.1 max_posV], v_ssi(vok), c_ssi(vok)); % fit sorted voltage
   
   % double boltzman for steady-state current
   inact.ss2_imax=fs2(2); inact.ss2_k=fs2(3); inact.ss2_v05=fs2(4);
   inact.ss1_imax=fs2(5); inact.ss1_k=fs2(6); inact.ss1_v05=fs2(7);
   inact.ss2_err=sqrt(sum((fit_ssi2-c_ssi(vok)').^2))/length(fit_ssi2); inact.ss2_Er=-81.5; inact.ss2_t=T;
   % single boltzman for peak current
   inact.ss0_imax=fs(2); inact.ss0_k=fs(3); inact.ss0_v05=fs(4);
   inact.ss0_err=sqrt(sum((fit_ssi-c_ssi(vok)').^2))/length(fit_ssi); inact.ss0_Er=-81.5; inact.ss0_t=T;
   inact.ss_V=v_ssi; inact.ss_I=c_ssi;
   
   watchoff;
   fit_t=toc;
   QueMessage(sprintf('Inactivation IV calculations took %8.3f sec', fit_t));
   CONTROL(sf).inact=inact; % store inactivation results in its own structure
end
%ctl=CONTROL; % return the modified structure.

if(tailanalysis)
    QueMessage('Tail current analysis proceding');
    tail = tail_fit(ts1+ts2, ts3, 3); % anzlyze the tail currents...
    CONTROL(sf).tail = tail;
    
end;

if(plot_flag < 0) % in case desire NO plots
   return;
end;


% now we can create the figure

% defaults
downsample=1;
fsize=7;
msg_pos = [0.37 0.00 0.15 0.07];

% window control
h = findobj('Tag', 'VC_IV'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'VC_IV', ...
      'MenuBar', 'none', ...
      'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
      'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
      'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
      'NumberTitle', 'off');
   
end

figure(h); % otherwise, select it
clf; % always clear the window...

%Command Menu
uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''VC_IV''));' );
uimenu('Label', '&Print', 'Callback', 'print;');
uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
uimenu('Label', '&Export Analysis', 'Callback', 'export(''act'', ''pk_V'',  ''pk_I'');');
set(h, 'Name', sprintf('Voltage Clamp Analysis - File: %s', CONTROL(sf).filename))
orient landscape

%plot currents
subplot('Position',[0.07,0.40,0.42,0.50])
[recs,pts]=size(current_s);
%make time array
plot(time(1,:)',current_s');
set(gca, 'FontSize', fsize);
ylabel('nA')
u=get(gca, 'YLim');

%mark peak times
if(bitget(analysis, 1))
   tp=p1A:1:p2A;
else
   tp=p1I:1:p2I;
end

yp=tp*0+u(2);
line(tp,yp,'LineWidth',3);
% mark the average times
if(bitget(analysis, 1))
   tm=m1A:1:m2A;
else
   tm=m1I:1:m2I;
end

ym=tm*0+u(1);
line(tm, ym, 'LineWidth', 2);

set(gca,'XTick',[ ])
title(sprintf('%s   R[%i - %i]',DFILE.filename,DFILE.frec,DFILE.lrec));
v=get(gca, 'XLim');
v(2)=tmax;
set(gca, 'Xlim', v);
datac_setcrosshair(gca, 'VCIV_info', 'ms', 'nA', msg_pos);

%plot voltage commands
subplot('Position',[0.07,0.2,0.42,0.12]);
[recs,pts]=size(VOLTAGE);
plot(time(1,:)',VOLTAGE');
set(gca, 'FontSize', fsize);
ylabel('mV');
xlabel('ms');
set(gca, 'Xlim', v);
datac_setcrosshair(gca, 'VCIV_info', 'ms', 'mV', msg_pos);

%plot analysis results
subplot('Position',[0.58,0.5,0.4,0.4]);
hold on;
mi=0; ma=0; vi=0; vm=0;
if(bitget(analysis,1))
   plot(volt_ss_A,curr_ss_A,'ro','MarkerSize', fsize*0.7);
   plot(v_ssa, fit_ssa, 'r');
   plot(volt_trans_A,curr_peak_A,'m+', 'MarkerSize', fsize*0.7);
   plot(v_tra, c_tra,'gx', 'MarkerSize', fsize*0.7);
   plot(v_tra, fit_pka,'g');
   
   mi=min([min(curr_peak_A'),min(curr_ss_A'),min(c_tra')])*1.05;
   ma=max(curr_peak_A)*1.05;
   vi=min(volt_trans_A)*1.1;
   vm=max(volt_trans_A)*1.1;
    if(vi == vm) 
        vm = vi + 1;
    end;
    if(mi == ma)
        ma = mi + 1;
    end;
    axis([vi, vm, mi, ma]);
   ylabel('nA');
   ylabel('nA');
   [h,hobj]=legend('mean','boltz','peak','peak-mean','boltz', 0);
   %   set(hobj(1), 'Fontsize', fsize);
   set(gca,'Fontsize',fsize);
   datac_setcrosshair(gca, 'VCIV_info', 'mV', 'nA', msg_pos);
end
if(bitget(analysis,2))
   plot(volt_ss_I,curr_peak_I,'kx', v_ssi(vok), fit_pki, 'k', v_ssi(vok), ...
      fit_pki2, 'r', v_ssi(vok), c_ssi(vok), 'bo', v_ssi(vok), fit_ssi, 'b');
   ylabel('nA');
   mi=min([min(curr_peak_I'),min(curr_ss_I'), mi])*1.05;
   ma=max([curr_peak_I', ma])*1.05;
   vi=min([volt_ss_I', vi])*1.1;
   vm=max([volt_ss_I', vm])*1.1;
   axis([vi, vm, mi, ma]);
   ylabel('nA');
   %h=legend('mean','peak','peak-mean', 0);
   set(gca,'Fontsize',fsize);
   datac_setcrosshair(gca, 'VCIV_info', 'mV', 'nA', msg_pos);
   
end
hold off;
grid

%plot text stuff
subplot('Position',[0.55,0.05,0.45,0.4])
r_uncomp = (1-number_arg(CONTROL(sf).comp/100))*number_arg(CONTROL(sf).access);
axis([0,1,0,1])
axis('off')
nli = 12;
pfsize=8;
text(0,txtpos(1, nli),sprintf('%-12s R[%d:%d]   %-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
text(0,txtpos(2, nli),sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', pfsize);
text(0,txtpos(3, nli),sprintf('Ra %4.1f MO Comp %3.0f%% => Rc=%4.1f MO',number_arg(CONTROL(sf).access), 100*number_arg(CONTROL(sf).comp),r_uncomp), 'FontSize', pfsize);
text(0,txtpos(4, nli),sprintf('Ileak: %6.1f nA (Vhold=%6.1f mV)',CONTROL(sf).iHold,CONTROL(sf).Vm0), 'FontSize', pfsize);
if(bitget(analysis, 1))
   text(0,txtpos(5, nli),sprintf('Mean@ %4.1f-%4.1f ms,Peak@ %4.1f-%4.1f ms',m1A,m2A,p1A,p2A),'FontSize',pfsize);
   text(0,txtpos(6, nli),sprintf('Mean Boltzmann: z=%7.2f (k=%7.3f) V_{0.5}=%6.1f mV, E=%8.2f', ...
      act.ss_k, 1/(F_RT(T)*act.ss_k), act.ss_v05, act.ss_err), 'FontSize', pfsize);
   text(0,txtpos(7, nli),sprintf('Peak Boltzmann: z=%7.2f (k=%7.3f) V_{0.5}=%6.1f mV, E=%8.2f', ...
      act.pk_k, 1/(F_RT(T)*act.pk_k), act.pk_v05, act.pk_err), 'FontSize', pfsize);
end
if(bitget(analysis, 2))
   text(0,txtpos(8, nli),sprintf('Mean@ %4.1f-%4.1f ms,Peak@ %4.1f-%4.1f ms',m1I,m2I,p1I,p2I),'FontSize',pfsize);
   text(0,txtpos(9, nli),sprintf('SS Boltz (single): z=%7.2f (k=%7.3f) V_{0.5}=%6.1f mV, err=%8.2f', ...
      inact.ss0_k, 1/(F_RT(T)*inact.ss0_k), inact.ss0_v05, inact.ss0_err), 'FontSize', pfsize);
   text(0,txtpos(10, nli),sprintf('PK Boltz (single): z=%7.2f (k=%7.3f) V_{0.5}=%6.1f mV, err=%8.2f', ...
      inact.pk0_k, 1/(F_RT(T)*inact.pk0_k), inact.pk0_v05, inact.pk0_err), 'FontSize', pfsize);
   text(0,txtpos(11, nli),sprintf('PK Boltz (dbl, 1): z1=%7.2f (k=%7.3f) V_{0.5}=%7.1f mV, err=%8.2f', ...
      inact.pk1_k, 1/(F_RT(T)*inact.pk1_k), inact.pk1_v05, inact.pk2_err), 'FontSize', pfsize);
   text(0,txtpos(12, nli),sprintf('PK Boltz (dbl, 2): z1=%7.2f (k=%7.3f) V_{0.5}=%7.1f mV', ...
      inact.pk2_k, 1/(F_RT(T)*inact.pk2_k), inact.pk2_v05), 'FontSize', pfsize);
end

text(0.75,0.08,sprintf('%s',date),'FontSize',6);

orient landscape
drawnow
if (plot_flag > 0)
   print
end
if (plot_flag == 2)
   close
end

function [yp] = txtpos(lineno, nlines)

spg = 1/(nlines+1); % spacing
yp = 1-lineno*spg; % relative y position
return;
